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Random unitary matrices find a number of applications in quantum information science, and are 
central to the recently defined boson sampling algorithm for photons in linear optics. We describe 
an operationally simple method to directly implement Haar random unitary matrices in optical 
circuits, with no requirement for prior or explicit matrix calculations. Our physically-motivated and 
compact representation directly maps independent probability density functions for parameters in 
Haar random unitary matrices, to optical circuit components. We go on to extend the results to 
the case of random unitaries for qubits. 


The development of the boson sampling problem m 
has motivated fresh interest in studying Haar random 
unitary matrices (HRUs) [5] realised with optical circuits 
to act on multiphoton states. Simultaneously, develop¬ 
ments in integrated optics [THU] now facilitate the con¬ 
struction of large-scale optical circuits capable of actively 
realising any unitary operator |16j including HRUs. Fur¬ 
thermore, HRUs play an important role in various tasks 
for quantum cryptography HU and quantum information 
protocols [M [IS], as well as the construction of algo¬ 
rithms m- 

Here we present a simple procedure for choosing a HRU 
on an optical circuit, implemented in terms of recursive 
decompositions of a unitary operator [SUES], by choos¬ 
ing values of the physical parameters independently from 
simple distributions. This procedure is useful for appli¬ 
cations where the exact unitary description of the imple¬ 
mented circuit is less important than a guarantee that 
it is drawn from the correct distribution. While simi¬ 
lar parameterisations exist in the mathematical literature 
[23j , an operational application within linear optics is not 
widely appreciated. We extend the result to systems of 
qubits, by deriving a mapping between a linear-optical 
circuit on TO = 2" modes and a circuit operating on n 
qubits. Note that constructions for pseudo-HRUs on qu- 
dit and qubit systems are also available, serving as a 
general framework to investigate randomising operations 
in complex quantum many-body systems |24[ I25j . 

Choosing a HRU is analogous to choosing a random num¬ 
ber from a uniform distribution, in that it should be un¬ 
biased. The probability of selecting a particular unitary 
matrix from some region in the space of all unitary ma¬ 
trices should be in direct proportion to the volume of 
the region as defined by the Haar measure, which is the 
unique translation-invariant measure on the space of uni¬ 
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FIG. 1: Recursive decomposition of a unitary in the trian¬ 
gular scheme, (a) An m x m unitary transformation can be 
factored as a product of m unitary transformations Ri, each 
acting on a successively larger subspace (the subscripts 0, ..., 
m — 1 on the left label the modes of the transformation), (b) 
A linear optical Rm can be constructed from a cascade of 
beamsplitters and phase shifts on optical modes. Both the 
Cartesian basis, x and the physical basis, r are illustrated. 


tary matrices. As argued in reference |26j . the columns 
of an TO-dimensional HRU may be made up from vectors 
= {i'll "(^ 2 , ■ ■ ■ ) I'm} that are successively drawn from 
the unbiased distribution of unit vectors in the subspace 
of (to — * J- 1) dimensions, orthogonal to all previous vec¬ 
tors. The problem of choosing HRUs thus reduces to the 
problem of recursively choosing such a set of orthogonal 
vectors. 
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FIG. 2: Direct dialling of a HRU in a triangular linear optical circuit, (a) A 6 x 6 unitary operator implemented with a 
triangular six mode linear optical circuit, (b) The pdfs from which beamsplitter reflectivities should be chosen to directly 
implement a HRU. Those on higher rows are chosen according to polynomials with increasing bias towards lower reflectivities, 
(c) A variable reflectivity beamsplitter can be effectively implemented with an MZI composed of a phase shift, 9 between two 
1/2 reflectivity beamsplitters. The MZI phases may be chosen directly from the distributions shown, with increasing bias 
towards 9 = n. The integrated optics implementation with directional couplers results in a bias towards 9 = 0. Line colours in 
(b) and (c) correspond to shading in (a). 


As we will show, this approach is particularly relevant 
to recursive circuit decompositions that allow any uni¬ 
tary matrix to be implemented over m optical modes, by 
choosing appropriate values for beamsplitter reflectivi¬ 
ties and phase shifters. We hrst consider the triangular 
scheme US] shown in figures [H which is a variant of 
that proposed by Reck et al. [21]. and which represents 
an TO X TO unitary matrix U as a product of unitary op¬ 
erators labeled R„, U = El]- Each block 

R„ is chosen to transform the mode j = to — n, or the 
corresponding basis state, denoted by \'$rn-n), into the 
n-dimensional unit vector |fn), over modes j = m — n to 
k = m — 1, i.e., 

I^n) = Rn l^m—n) ■ (1) 

This vector undergoes further transformations under sub¬ 
sequent blocks Ri (n < i < to) to finally produce |/„) 
that occupies all to modes: 

l/n) = Rm • • -Rri+l |fn) • (2) 

Orthogonality between each of the m\fi) vectors is guar¬ 
anteed. Further, if the vector \vn) is chosen from the un¬ 
biased distribution of unit vectors in n dimensions, the 
property of left invariance ensures that |/„) does not be¬ 
come biased by the operation of the subsequent R^. 

The next and main task is therefore determining how 


an unbiased vector in n dimensions may be implemented 
with R„ by choosing values for the linear optical com¬ 
ponents from which it is constructed, according to the 
expansion shown in hgurej^b). To achieve this, consider 
the complex Gaussian vector in n dimensions: 


n—1 n—1 

|v„) = '^Zi IT-J = ^ , (3) 

i—Q i=0 


where IT^) denotes the ith basis state and the Zi are 
independent and identically distributed normal random 
variables with the probability density function (pdf), 

Vziiz) = I/tt exp . This independence means 

that the pdf for v„ is the product of the pdfs for these 
elements and depends only on the magnitude of the vec¬ 
tor: 


(x) = J_e-^^o+^i+-+=Or.-i) ^ ^g-|v.p 

!■ ^ z _ T). _n. ’ 


( 4 ) 


where Xi = \zi\^. We now show how this basis x, which 
we call the Cartesian basis, can be mapped to a new 
basis, r. We call the latter the physical basis, since, as 
we demonstrate below, it contains the variables corre¬ 
sponding directly to components in a physical realisation 
of the vector in linear optics. Namely, we denote by tq 
the power of the input to the given block R„, while the 
other ri stand for the reflectivities of beamsplitters (see 
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also figure [T|b)). Next, combining the definition Q and 
equation we find 

zo = (5) 

i 

Zi = y/ror,+i \/l - ru 0 <i<n-l, (6) 
k=l 

where the matrix (in the Pauli basis) B{r) = ^pra^ + 
yr — ra^ has been used to describe a beamsplitter as a 
function of its reflectivity. Finally, taking into account 
that Xi = I Zip, we find, 

(7) 

0 < i < n — 1 (8) 

(9) 

We must show that the pdfs for the vector v„ are separa¬ 
ble in the physical basis so that the experimental parame¬ 
ters can be chosen independently. We also need to derive 
the form of the marginal distributions for the ri and tpi, 
from which experimental parameters must be chosen to 
obtain a Haar unitary. Since there is no functional depen¬ 
dence on the at parameters in equation (|^ and there is a 
one-to-one mapping ai ^ p, these phases can be chosen 
uniformly and independently from the interval [0, 27r). 

Finding the pdfs for the beamsplitter reflectivities re¬ 
quires a more careful change in bases, using the Jacobian, 

= 7^v„(x) |det J(x,r)|. (10) 

The pre-factor from Q is expressed in the r basis sim¬ 
ply as exp(—ro), so is trivially separable. We therefore 
consider the Jacobian matrix 

= ( 11 ) 

with 


n—1 

ro = y^^Xk 

k=0 

X^-l 

n = —- 

P,k=i-1 

■ 


where the variable Cn = 1 has been introduced for con¬ 
venience. 


We show that this form of matrix (lower Hessenberg) can 
always be transformed into a lower triangular matrix— 
for which the determinant is simply the product of the 
diagonal elements—by elementary operations, which do 
not change the absolute value of the determinant. 


The hrst step is to perform a set of operations on the 
j = 0 column, Cq, that set the upper n—1 terms to zero, 
as follows: 


T(fc— 1) 

( fc ) ( fc - 1 ) "’ fc - 1,0 

Cq =c(, '-Cfc--, 

Jfc—l.fe 


■"0 


(18) 


where k runs from 1 to rn — 1, and J^^g are those 
quantities after k operations, and is the fcth column. 
We can then place the column Cq as the rightmost col¬ 
umn, at which point the matrix is lower triangular. After 
the procedure is complete the element Jo,n-i = 1 (see ap¬ 
pendix for detailed proof). 


The Jacobian determinant is given by multiplying the 
diagonal elements of the shifted matrix 


n—1 

det J(x,r) = J'j^i (19) 

i=l 


which are given by equation (16). 
the pdf in the r basis is 


The explicit form of 


T^vjr) 


k^l 


( 20 ) 


which is manifestly separable. 


It can be verified by explicit integration that this expres¬ 
sion is appropriately normalised. Since the pdf is sepa¬ 
rable in this basis, the variables are independent, and 
can be chosen according to their marginal distributions. 


xo = ro ri ( 12 ) 

i 

Xi = ro r ^+1 JJ(1 - rfc) 0<i<n-l. (13) 

For the four cases 


I 


= Vi+I 11 (1 - Tk) 

j = 0 

(14) 

k^l 



■1 k=l 

0 < j < 1 

(15) 

=^0 n (^ “ ’’fc) 

j = i + l 

(16) 

k^l 



= 0 

j>i + l 

(17) 


Vr^Pr) = {n — i) {I — rp * ^ 1 < i < n, (21) 

where, for clarity, denotes the reflectivity of the zth 
beamsplitter in the nth rotation, R„. We now inte¬ 
grate over ro to obtain a compact form for the pdf of 
n-dimensional unit vectors, 

n—1 

(r) = (n - 1)! (1 - r„,fe)"-'=-^ (22) 

fc=i 

and express the pdf for the full circuit of beamsplitters, 
Vc{x) as the product of the pdfs for the diagonal arrays 
of beamsplitters: 


rc(r) = n 


1-1 


(■? “ ^)'n (1 “’’ja) 


j-k-l 


fe=l 


(23) 
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Recalling the beamsplitter transformation B{r) = 
V^o-z + Vl — r(7x, we note that a variable reflectivity 
beamsplitter can be constructed as a Mach-Zehnder in¬ 
terferometer (MZI), from a variable phase shifter 9 be¬ 
tween two 1/2 reflectivity beamsplitters, H = 5(1/2), to 
give By(9) = cos |/ -|- i sin |era, (up to a global phase). 
It is then useful to re-express the pdfs in terms of MZI 
phase shifts. The further change of variables, r = cos^ |, 
gives 



VfM = in 




2{n—i) — l 


(24) 


In the setting of integrated optics, where beamsplitters 
are implemented with directional couplers on waveguides 
according to 5(1/2) = -^(I + iax) for reflectivity of 1/2, 
the pdfs are given by ( |24[ ) but with sin and cos functions 
interchanged, i.e.. 




9 

COS - 
2 


2{n—i) — l 


(25) 


FIG. 3: A 6 X 6 unitary operator implemented with a 

six-mode linear optical circuit according to the rectangular 
scheme. Here, f„,i stands for the reflectivity of the ith beam¬ 
splitter of the block Within each R„ {n = 2,..., 6), 

we enumerate the beamsplitters according to the sequence 
s, which consists of n — 1 indices, with odd (even) numbers 
arranged in descending order and followed by even (odd) num¬ 
bers, arranged in ascending order (see also main text). In this 
figure, the beamsplitters in the ith row mix the modes i — 1 
and i (e.g., the beamsplitters of the third row, f 4 , 3 , fep and 
f 5 , 2 , couple the modes 2 and 3). 


In practical terms, an optical circuit composed of beam¬ 
splitters and variable phase shifters can directly dial up a 
configuration corresponding to a HRU, by choosing phase 
shifter values from the derived pdfs. A six mode example 
is given in figure 


We note that the version of the triangular scheme used 
here, in which each beamsplitter in every block R„ 
couples two adjacent modes, differs from the original 
scheme m, in which the first mode is consecutively 
coupled with modes 2,3, ...,n. It is easy to check, 
however, that the mapping 0-0 can be applied to 
the original scheme as well, by replacing r with 1 — r 
and relabelling the output modes: {xo,Xi, > 

{xm-i^xo ,..., Xm- 2 }- Such a change of variables does not 
affect the Jacobian determinant and the final expression 
for the reflectivity pdfs for the original scheme is obtained 
by replacing r with 1 — r in equation (21) (the phases are 
again chosen uniformly and independently from the in¬ 
terval [0, 27r)). 


Next, we analyse the alternative decomposition of uni¬ 
tary matrices, proposed by Clements et al. [52], which 
corresponds to a rectangular mesh of beamsplitters and 
phase shifters, as shown in figure]^ for six modes. While 
the triangular scheme might be more resilient to loss and 
other errors in experiments in which only a small pro¬ 
portion of its (upper) input ports are accessed, the rect¬ 
angular scheme is likely to be beneficial for experiments 
that involve accessing most of its inputs. The more com¬ 
pact rectangular scheme may also fit a greater number 
of modes on standard wafers used in the fabrication of 
integrated photonic circuits. 


The rectangular scheme obeys the blocked structure, 
analogous to the triangular scheme described above. 
That is, an TO X TO unitary matrix U can be written down 
as a product of blocks R„ (hereafter the tilde refers to the 
decomposition of reference 1221). Each of these blocks, as 
previously, transforms the mode m — n into a vector over 
modes to —n up to to — 1 (see also figure 0b)). More pre¬ 
cisely, for odd (even) to, U = llY=i ^^j-i Rm- 2 i 

(U = Moreover, the 

mapping of equations Q -0 for the operator R„ for the 
triangular scheme can be used for R„ as well, by a simple 
substitution. Namely, for even (odd) to we replace fn,i 
by rn,s(i), Vi, where s is a sequence of n — 1 indices, with 
odd (even) numbers arranged in descending order and 
followed by even (odd) numbers, arranged in ascending 
order (e.g., for to = n = 6, we have s = {5,3,1, 2,4}). 


This substitution leaves the corresponding Jacobian de¬ 
terminant unaffected. Therefore, the pdfs given in equa¬ 
tion (21) for reflectivities for the triangular scheme 


correspond to that of the rectangular scheme, but for 
In other words, ^{r) = T’f^^(.j(r). Subse¬ 
quently, we find 


'Pfy.iir) = [n- s(i)](l - f) 


n—sfi) —1 


(26) 


Alternatively, one can reorder the reflectivities fn,i ac¬ 
cording to the sequence s, as is done in figure]^ yielding 
Vry i (r) = Vfy^i {r). Finally, the phases of the rectangular 
scheme, analogous to the triangular scheme, are chosen 
uniformly and independently from the interval [0, 27r). 


Given the above parameterisation of Haar-random opti¬ 
cal circuits, we now address the effects of errors, caused 
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FIG. 4: Coverage cov^ of the unitary space versus the circuit 
size m. The phase shifters are assumed to cover their full 
range [0,27r), while the range of reflectivities is restricted to 
[|e|,l — |e|], where random errors e are drawn from a zero- 
mean normal distribution. The curves correspond to different 
variances a of the errors (a = {1,5,10,20} x 10“^ from the 
upper to the bottom curve). For each m, the coverage is 
averaged over many realizations of e. 


by imperfections in integrated photonics manufacturing. 
Before going into detail, we emphasize the important fea¬ 
ture of our approach: due to the separability of the de¬ 
rived probability distributions, errors on a given compo¬ 
nent of the circuit do not propagate to other indepen¬ 
dently chosen parameters. In turn, a major source of in¬ 
dividual errors is the imperfection of directional couplers. 
Used to implement the balanced beamsplitters of MZIs, 
directional couplers should ideally couple 1/2 of the light 
between waveguides so that each MZI can achieve the 
full reflectivity range. Fabrication tolerances, however, 
introduce errors and limitations on this range. Further¬ 
more, we note that upper MZIs in the triangular scheme 
and central MZIs in the rectangular scheme are those 
most sensitive to errors, according to their polynomially 
growing pdfs. 


Although schemes exist to minimise the effect of such 
errors and produce near perfect MZIs HZIHH] it is worth¬ 
while considering the influence of many small errors over 
a large circuit (this simple model is also useful to the 
qubit picture that we develop below). As an estimate 
to this effect we address the range of unitary operations 
covered by the proposed parameterisation, which we eval¬ 
uate in terms of the coverage of the unitary space (see 
references [ 23121 ] for more details), 


COVm 




'Id 


/u^U 


(27) 


That is, cavrn is the ratio between the reachable and full 
unitary spaces, assuming that the phase shifters cover 
their full range [0,27r). The range of MZI reflectivities, 


in turn, is [|e|, I — jej], where e is a small random error. 
In figure]^ we plot the coverage versus the circuit size to, 
which shows that for such moderate errors our parame¬ 
terisation achieves high coverage rates. Since the pdfs for 
the triangular and rectangular schemes have been shown 
to be equivalent and independent, the coverage plotted 
in figure]^ is valid for both. 

We now briefly show how these results may be ex¬ 
tended to the scenario of quantum information process¬ 
ing with qubits, independently of any particular phys¬ 
ical implementation. We suggest a mapping between 
a unitary operation on to = 2^ optical modes and 
the same unitary operation on p qubits, such that the 
pdfs derived above can be directly applied to systems 
of qubits. Labelling the optical modes as qubit basis 
states {|0...00), |0...0I), |0...10),..., |I...1I)} we map the 
optical beamsplitters and phase shifters to single qubit 
Hadamard gates, and n-qubit logic gates where the state 
of a single target qubit is transformed depending on the 
states of n — I control qubits. The target qubit opera¬ 
tions are the NOT gate or qubit-flip operator, a^, and 
the qubit-phase gates, $ = and $ = Each 

optical phase is mapped to a n-qubit $ or $ logic gate, 
and each optical beamsplitter is mapped as an MZI to 
a n-qubit $ or $ logic gate between two single qubit 
Hadamard gates on the respective target qubit. 

The mapping can be understood with reference to fig¬ 
ure which explicitly details the case for 3 qubits and 
8 optical modes and present a full circuit example for 2 
qubits and 4 optical modes. The target for the n-qubit 
phase operations is always the final qubit; the condition¬ 
ing configuration of the control qubits determines which 
element in the qubit space receives the phase. The addi¬ 
tion of Hadamard operations on the final qubit allows the 
mapping of 1/2 reflectivity beamsplitters, and therefore 
MZIs, that operate between pairs of optical modes that 
differ in labelling only by the final bit. The further ad¬ 
dition of n-qubit NOT gates allows MZIs to be mapped 
from pairs of optical modes that may differ in labelling 
by more than one bit. Any subset of the MZI opera¬ 
tions may be implemented on qubits by simply omitting 
controlled phases where appropriate. 

While not designed to be optimal, this one-to-one map¬ 
ping between n-qubit phase gates and optical MZIs il¬ 
lustrates one way in which the distributions expressed in 
figure [^c) may be used to directly implement a HRU on 
qubits. 

We have presented a recipe to directly generate HRUs in 
linear optics with a proof that is straightforward in com¬ 
parison to previous works Iii23|. Experimental confor¬ 
mation of these results can make use of tomography that 
does not require further optical circuitry [30]. The for- 
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(a) 




FIG. 5: Mapping a linear optical circuit to a unitary on qubits. An empty or solid circle indicates that the operation is 
conditional on the ‘0’ or ‘1’ state of the control qubit, respectively, (a) A 2-qubit and 4 optical mode example, (b-d) The 
elementary operations for a 3 qubit (8 optical mode) unitary, (b) Phase shifts on single elements in the qubit space result from 
a judicious choice of conditions for the control qubits, (c) The addition of Hadamard gates allows beamsplitter operations to 
be implemented on elements in the qubit space that differ only in the state of the final qubit, (d) The further addition of 
n-qubit-NOT gates allows the mapping of beamsplitters for elements that may differ in the state of more than one qubit. 


mula in its general form is applicable to boson sampling 
where Haar unitaries are required, and the extension to 
systems of qubits invites wider applications. 
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one of the authors, Nick Russell. Those of us who knew 
him are grateful for his contributions to our work and 
our lives, which we continue to miss. We acknowledge 
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search Council (EPSRC), the European Research Coun¬ 
cil (ERC), including QUCHIP (H2020-FETPROACT-3- 
2014: Quantum simulation), the U.S. Army Research Of¬ 
fice (ARO) grant W911NF-14-1-0133. A.L. acknowledges 
support from an EPSRC early career fellowship. 
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Appendix: Converting a Lower Hessenberg matrix 
to a lower triangular matrix 

We set Q = 0 with column operations by subtracting 
column Cfc multiplied by an appropriate scalar: 


c 


(fe) 

0 


= c 


(fe-1) 

0 




t(^—1) 

J/c—l./c 


(A.l) 


The effect on all the other elements of Cq is to remove 
the dependence on r*,, which we can prove inductively. 

Suppose that after k such operations, the upper k el¬ 
ements of Cq have been set to zero and the remaining 
elements have no dependence on r; for 0 < I < k. We 
can express the elements of as: 


Rfe) _ / r,+i riLfc-m (1 - n) , i>k 
1 0 , i<k 


The base case is k = 0, where the expression in (14) 
corresponds to this general form. We now perform the 
{k -b 1) th operation on all non-zero rows (i.e. i > k 1): 


t(^) 

-[(fc+l) _ -,{k) '’k,0 -[ 

"^i.O ~ '^i.O T 'Jj.fc+l 

i 

= r^+l {l-ri) + 

l^k-\-l 
k 

rk-\-i n 

l=k+l roTi+i -r-r , 

li(l-G) 


ro 




(1 ’’fc-l-l) 


1^1 


= r. 


4 + 1 


n 


'I’fc+l'l’i+l 


n (i-’’*) 


Z=fc+1 ^ ^ ^ /=fe + l 

i 

= r,+i (1 - rfe+i) (1 - n) + 


/= fe +2 




n 

l^k+2 

i 

= Tj+i (1 - Tfe+I-b rfe+i) (1-n) 


Z=fc+2 


> i+1 


n 


Z = fc-|-2 


We recover the expression in (A.2), thus proving the re¬ 
sult. After n — 1 iterations, we find that o = = 1) 

recalling that r„ = 1 was a variable introduced for con¬ 
venience. 








